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a projectile, an ace-old problem occurs. A conflict between 
optical quality and aerodynamic drag exists. A new approach 
for solving this problem is being developed using gradient 
refractive index, GRIN, for lens design. The equations for 
lens design incorporating GRIN are formulated for the special 
case of spherically symmetrical GRIN. A computer code for 
designing a GRIN lens and for determining aberrations has 
been developed. A different program (LENS) calculates lens 
design and lens performance for the case of homogeneous 
refractive index. Results of this program provide a check 
on the more complex program which includes GRIN. The 
computed program, calculates the spot diagram, for both mer- 
idianal and skew rays. The spot diagram, which is the inter- 
section of rays with the focal plane, is plotted, and the a r 
and centroid location of the spot diagram are calculated. 

Also a cumulative energy diagram is obtained from the program. 
The computer code has been organized to calculate the items 
enumerated above for either a circular cone for the inside or 
outside surface of the pointed lens. For homogeneous 
refractive index, both cases are included in the program 
discussed in this report. However, for the GRIN case, only 
the complete program for a lens with a conical surface on 
the outside is reported here. The alternate case is being 
developed by a thesis student. 
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Features of the computer program include the following: 
use of many subroutines for clarity and for easy modification, 
and use of program structure for easy integration into an 
optimization code. 

Sample lens designs and lens performance are presented. 

The computer code provides the tools for lens design. Consider- 
able additional work, is needed for selecting the best lens 
design. 
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1. INTRODUCTION 



A long-term conflict exists between optics and aero- 
dynamics for missile nose design. Optics requires a blunt 
shaped nose similar to that used on the Sidewinder missile. 
Favorable aerodynamics for supersonic missiles requires a sharp 
nose which permits an attached shock wave. For rocket 
powered missiles or projectiles, the sensor and associated 
optics occupies the complete nose area. The penalty of a blunt 
nose is high drag. For ramjet propelled missiles or pro- 
jectiles, the penalty of a blunt nose is poor pressure recovery 
by the inlet; this statement assumes a nose inlet is used. 

In times past, the lens designer had available three design 
parameters; a satisfactory lens in the general shape of a cone 
has not been designed. A lens can be designed using homo- 
geneous refractive index which focuses rays; however, off-axis 
rays have very large aberrations. The thesis by Terrell [1] 
discusses the design of a conical lens. 

A fourth design parameter has become available which is 
spatially variable refractive index. The associated technology 
is known as GRIN; see references [2-13]. GRIN is short notation 
for Gradient Refractive Index. The thesis by Frazier [14] 
briefly introduced the use of GRIN in lens design; likewise 
the thesis of Terrell [1] has a brief discussion. The thesis 
of Carr [15] will have more extensive results. 

The development of materials which have variable refractive 
index has been motivated by optical fiber technology. Optical 
fibers can propagate a single mode with low losses if an ap- 
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propriate spatial variation of refractive index exists. Based 
on the information in reference [2] , the maximum variation 
of GRIN is limited to a few percent. The report contains 
twelve sections. In section 3 the problem of the homogeneous 
case is solved (program LENS) while the GRIN lens design is 
formulated in section 4 (program GRIN) . The flow charts to 
the conputer programs are given in section 5 and 6, and a 
user's guide to each of the programs is illustrated in section 
7. The programs are listed in sections 9 and 10. 

2. BRIEF OVERVIEW OF PROBLEM FORMULATION 
Define the axial distance as x and radial distance as 
r in a cylindrical coordinate system. The lens is to be 
axisymmetric so that angle is not a variable. 

In general the shape of the front surface, is a function 

of r 

x = F f (r) (2.1) 

where subscript f denotes front surface. Likewise the rear 
surface is a function of r 

x = F r (r) (2.2) 

where subscript r indicates the rear surface. 

As a simplification, either the front or rear surface is 
a circular cone. Terrell [1] treated the case of a conical 
front surface; Carr [15] is investigating the case of a conical 
rear surface. This report develops a computer code which 
calculates either case. A single consistent computer program 
is vised for both cases; this fact is important for application 



2 



of optimization codes. 

In general, the GRIN can be specified as 

n = F (x, r) (2.3) 

The function F involves both spatial variables. However, a 
less complex GRIN is used which is spherically symmetrical 
refractive index. In this case 

n = n (x c ,r' ) (2.4) 

where r' is radial distance from the center of symmetry for 
the sphere. The location of the center of the spherical GRIN 
is on the lens axis at axial position x c . 

One very important computational advantage of spherical 
GRIN is that the rays within the lens are confined to a single 
plane. In the case of skew rays, the plane of the rays is 
not a meridianal plane of the lens. 

The equation for a ray trajectory in a spherically 
symmetric GRIN geometry is 



0 = 
e 



dr ' 



/ ,2 2 



n ( r ' ) -3‘ 



(2.5) 



where E is a constant for a ray. The value of E is obtained 
from 

E = en Q (r^) r^sin^Q (2.6) 

Define r' as a unit vector in the direction of a line extending 
from x c to point P in Figure 1. Define R as a unit vector of 
a ray which has been refracted by the surface at point P. The 
angle Fg 



is obtained from 



( 2 . 7 ) 



cos 




•I 




CENTER OF SYMMETRY OF GRIN 



Figure 1. Geometry for Ray Tracing in a Spherically Symmetrical 
GRIN medium. The ray vector, R , and the radius 
vector, r' , are illustrated. 
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In equation (2.6), the e is a sign function having a value of 
+ 1.0 or - 1.0 depending on the value of f within the lens. 



(+ 1.0 ¥ < tt/2 

e =< " (2.8) 

(- 1.0 ¥ > tt/2 



Equation (2.5) can be integrated fcr certain functions n(r'). 
One function is 



n(r’) 



A + B(r'/r 0 ) 2 



h 



(2.9) 



where A and B are constants. Tq is a value associated with a 
specific lens design; a convenient value for r Q is that value 
for the ray defining the lens aperture. 

Combining equations (2.5) and (2.9) and completing the 
integration yields 



e (r’) = -35 
e 



sm 



-1 



? 9 

2E /r’ -A 
r 

(A 2 +4BE 2 /r Q V* 5 



>-sir. 



-1 



2 2 

2E Z /rV -A ' 
(A 2 +4BE 2 /r Q 2 ) 45 ; 



J (2. 10) 



Equation (2.10) is an explicit relation for 0 g as a function 

of r’ . r , r £ are the radii to the rear surface and to the 
r f 

front surface respectively. An explicit relation for r’ as a 
function of 9 g can be derived; the result is 



r '(r, f) (e e> " 



/2 | E 



A+ / A 2 +4BE 2 /r, 



sin / + 20 +sin 2 
e 



2 £Z/r (r, f)-aT 



A +4BE^/r 0 j 



( 2 . 11 ) 
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For details look at section 4.1. For certain phases of the 
calculation equation (2.11) is needed. 

Thinking ahead to numerical optimization, the function 
of equation (2.9) was chosen to avoid numerical integration 
of equation (2.5) . Equation (2.5) must be evaluated several 
times for each ray. A lens design may involve 100 to 1000 
rays. Numerical optimization of the lens requires calculation 
of many different lens shapes. Hence, evaluation of equation 
(2.5) may occur several thousand times in the process of lens 
optimization. The GRIN-related parameters to be determined 
in optimization are A, B and x c . 

3 . EQUATIONS F OR THE HOMOGENEOUS CASE 
3 . 1 Meridi a n Plane 

The homogeneous case is included here for several reasons. 
First, the homogeneous case is a simple example in which the 
intricacies of lens design are apparent in considerable clarity. 
Iteration is not required for design of a lens with a homo- 
geneous refractive index, whereas for GRIN, iteration is 
necessary. Second, the homogeneous case provides a meaningful 
check of the GRIN program. To conduct the check, the value of 
B is zero in equation (2.9). Obviously the answers from GRIN 
program with B equal to zero must agree with. answers from the 
homogeneous case. 

The book by Kingslake £ 16 J has been helpful in the program 
development . 
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Figure 3. Geometry and angle definitions for the homogeneous refractive index case 



3.1.1 The Eauations 

Refer to Figure 3 for definition of the symbols. 

n (J) + I 2 ' (J) =0 (I) (3. 1) 

sin (I 2 [ J] ) = n 32 • sin(I 2 '[J]) (3.2) 

The notation n^ is used for the ratio n 3 /n 2 * 

I 1 '(J) = U/2-ct) - (n (J) + l 1 (J)) (3.3) 

sin (I^[ J] ) = n 21 • sind^' [J]) (3.4) 

a + I 1 (J) + U = tt/2 (3.5) 

3.1.2 Solution 

The angles a, 8(J) are inputs and are known. 

From eq . (3.5): 



Ij_(J) = tt/2 - a - U (3.5.1) 

From eq. (3.4) : 

I 1 ' (J) = sin _1 [n 12 • sin(I 1 [J])] (3.4.1) 

Substitute eq. (3.5) to eq. (3.3): 

n (J) + I 2 (J) = I 1 (J) - I 1 ' (J) + U (3.3.1) 

But from eq. (3.2) : 

I 2 '(J) = sin * sin(I 2 [J])] (3.2.1) 

From eq . (3.3.1): 

I 2 (J) = [I 1 (J) - I 1 ’(J)] - ri(J) - U (3.3.2) 

Substitute eq. (3.3.2) into eq. (3.2.1): 

I 2 '(J) = sin -1 [n 23 • sin(l 1 (J) - I 1 ' (J) - n (J) + U) ] (3.2.2) 

Substitute eq. (3.2.2) into eq. (3.1): 

n (J) + sin -1 [n 23 • sin(l 1 (J) - 1 1 ' (J) - n (J) + U>] = 6 (J) 

(3.1.1) 



9 



Equation (3.1.1) is solved for n ( J) : 



n 23 • sin(I : (J) - I 1 ' (J) - n (J) + U) = sin[ 8(J) - 
Use the formula: 

sin (A ± B) = sinA cosB ± cosA sinB 
And substitute: 

x = sin(n(J)) ; = cos(n(J)) 

Equation (3.1.2) will turn into the form: 



n (J) ] 



x = 



v m? 



(3.1.2) 



(3.1.3) 



Where : 



A = 



n 23 ’ - 1^' (J) + U] - cos[B(J)] 



n 



23 



sin [I . ( J) - I ' (J) + U] - sin [ 3 ( J) ] 



(3.1.4) 



Therefore : 



n (J) = sin -1 [x] 



(3. 1.5) 



Knowing n (J) , I 2 ' (J) can be solved from equation (3.1) , and 
then, I 2 (J) can be solved from equation (3.2). 

3.1.3 Coordinates of the Front and of the Rear Surfaces 
3 . 1 . 3 . 1 The First Ray 

All the magnitudes will be normalized 
with respect to the radius R. R is the aperture for the lens. 
There fore : 



R 



1.0 



y ! ( i) 

X-J^l) 



R 

R 

tanToJ 



(3.6) 
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The thickness of the lens is defined: 



T ( 1) = J [x 2 (l) - x^l)] 2 + [y L ( 1) -y 2 (l)] 2 (3. 



Therefore, the appropriate coordinates of the rear surface are 



(3 



x 2 (l) = x 1 (l) + T ( 1) • cos[tt/2 - a - I^'(l)] 
y 2 ( 1) = y 1 (l) - T ( 1) • sin [tt/2 - a - 1^(1)] 

3 . 1 . 3 . 2 The Other Rays: Coordinates of the Rear 

Surface 

The new point of the rear surface 
x 2 (J),y 2 (J) will be the intersection between the surface and 
the ray. 



(a) The equation of the surface: 









y 2 (J) = m[x 2 (J) 


- BB' ] 




(3. 






m = 


tan[-r/2 - n ( J - 1) ] = 


: cot [n ( j 


- 1) ] 


(3. 








BB' = BF - B' 


F 






B’F = 


^2 


(J 


- 1) * [ cot [tt/2 - n ( J 


- D] + 


cot [B ( J - 1) ]] 




3 


1 F 


= y 


2 ( J - 1) • { tan [n ( J - 


1) ] + cot[S(J - 1) ]> 


(3. 


• 

• • 


^2 


(J) 


= cot [n ( J - 1) ] { x 2 ( J) 


- [BF - 


B ' F ] } 


(3. 


(b) 


The 


equation of the ray: 










^2 


(j) 


= -tan[B ( J) ]x 2 ( J) + 3F • tan 


[3 (J) ] 


(3. 






• 

• • 


y 2 (J) = tan[B(J)] * (- 


x 2 (J) + 


BF} 


(3. 



From equations (3.9.3) and (3.10.1), the coordinates of the 
next point at the rear surface are found; those are: 



v /-rv _ BF{tan [ 8 ( J) ] + cot[n (J - 1)]} - cotfn(J-l)] * 
“2 ° cotfnT j I) - ! ~~+ tan[S(J)l 

y 2 ( J) = tan[$(J) ]{-x 2 (J) + BF } 



7) 



8 ) 



9) 

9.1) 



9.2) 

9.3) 

10 ) 

10 . 1 ) 



B * F 

ID 
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3. 1.3. 3 Coordinates of the Front Surface 



(.a) Equation of the rear surface: 





y^(J) = (tana) * x^(J) 


(3.12) 




(b) Equation of the ray: 






y 1 (J) = m[x 1 (J) - x 2 (J)J + y 2 (J) 


(3.13) 


where : 


m = -tan[u/2 - a - I ' (J) ] 


(3.13.1) 


or: 


m = -cot [a + I. ' (J) ] 


(3. 13.2) 



From equations (3.12) and (3.13) the intercept point is 
cal culated: 

y 2 (J) + x 2 (J) • cot [a + I 1 ' (J) ] 

x, (J) = 

tana + cot[a + I.'(J)] 

(3.14) 



y 1 (J) = (tana) • x^J) 



3 . 1 . 3 . 4 The Thickness at the New Point 

T ( J) -J [x 2 (J) -x 1 (J)] 2 + [y 1 (J) - v 2 (J) ] 2 (3.15) 

3 . 2 Skew Rays in the Homogeneous Case 
3.2.1 Coordinate Transformation 

The x,y coordinates will be shifted by the 

angle a_ . 

P 

x' = x * cos (a ) - y • sin(a ) (3.16.1) 

y' = x * sin(a ) + y * cos (a ) (3.16.2) 

hr hr 

z' = z (3.16.3) 



3.2.2 Calculate the Intercept Point of the Front 
Surface 

Assuming that the (0,0) point is at point A 
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(Refer to figure 2.), the following equations can be derived: 

(3 

,2 



(a) Sa*: R ray 2 = Y 2 + z 2 



R 



ray 



y ' 

* 0 



cos (a T 
P 



- (x + AB) • tan (a ) 

P 



+ z, 



(3 



Where R is the radius of the ray and the values for y_ ' , 
ray ■‘0 

z 0 ' are known . 

(b) Cone: R _ = tan (a) * [x + AB] (3 

cone 

Solving equations (3.17.2), (3,13) for x results: 



XO = - _£2- S + Y. 



'CDS T + A 2 D 2 ♦ CZD 



A 2 D 2 



A 2 D 2 



(3 



YO = 



ZO = 



Y o’ 

— - - XO . tan (a ) 

cos (a ) p 

ir 

V 



(3 

(3 



Where : 



2 2 

A-D_ = (tana) - (tana ) 

2 2 p 

y 0 ' 

CDS = tana P + (AB)[A 2 D 2 ] 

p 



(3 

(3 



CZD = ( -£0- ) 
\ cosa / 



, 2 

y O \ . , ,2 



+ Z 



0 



2 2 V 

(AE) [A 0 D 2 ] - 2(AB) rr ^ r - tana. 



cosa 



3.2.3 Direction Cosines of the Normal 



(3 



f =y 2 + z 2 - (tana) 2 (x + ab) c = 0 



(3 



^ = -2 (tana) 2 ( x + AB) 



(3 



.17.1) 

.17.2) 

. 18) 

.19.1) 

.19.2) 

.19.3) 

.20.1) 
.20. 2) 

.20.3) 

. 21 ) 

. 21 . 1 ) 
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Figure 4. Diagram showing coordinate transformation for 
off-axis rays. 
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(3.21. 2) 



6f _ ~ . 6f _ 5_ 

~ ^ z 

6 V ^ 5 Z 



2,% 



| Af j = 2 [ (tana) 4 (x + AB) 2 + y 2 + z 2 ] 



Substituting equation (3.21) can turn equation (3.21.4) i 
the form: 

I “—►I 0 +- rv 

(x + SB) 



Ilf I = - 2tana 



cosa 

The direction cosines of the normal are therefore: 






I A f | 



= -sma 



Define : 



rv - (Sf/Sy) 

^ " lifl 



a = tan" 1 ( Z /Y) 



2Y 



2tana/cosa ( x + AB) 



y 2 + z 2 = (tana) 2 ( x + AB) 2 



y 2 (l + tan “a) = tan"a(x + AB) 2 



y = (cosa) (tana) (x + AB) 

The re f o re : 

= ( cosa) (cosa) 

Similiarly : 

n = l l i ^z .L = (sina) (cosa) 
Z I Af I 



3.2.4 Refraction at the Front Surface 

t 

The skew ray R will be defined: 

s = p x • <s x ) + ^ • (i. f ) + Rj • d z ) 



(3.21.4) 
to 

(3.21.5) 

(3.22) 

(3.23) 

(3.23.1) 

(3.23.2) 

(3.23.3) 

(3.24) 

(3.25) 
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Figure 5. Cross section of lens in a plane normal to the lens 
axis. Point P is the point where an off-axis 
ray intersects the front surface of the lens. 
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, R z are the direction cosines of the ray, 
and they have the value: 

R x = cosotp ; Ry = -sina p ; R z = 0 (3.26) 

The surface normal is defined by the following expression: 

A /V /v 

N = n x • (e x ) + n y . (e y ) + n z (e z ) (3.27) 

Where the direction cosines of the normal surface were 
defined previously (equations 3.22, 3.23.3, 3.24). 

The refracted angle at the front surface is defined: 



Whe re R x , 



e 



-l 

cos 



R . N 

Ir| . |n| 



(3.28) 



Whe re : 






N 





+ 




+ 




1 



(3.29) 



There fore : 

9 = ccs -1 {R x . n x + R,^ . n y + R 2 . n z ) (3.30) 

<p = tt - 9 

<j>' = sin 

3.2.5 Refraction Equation s 

Referring to figure 6 the refraction at the 
rear surface with regard to the front surface will be as 
follows [16, p. 147]: 



(3.31) 



n z R' = n^R + [n 2 Cos(<^: , ) - n-,cos(4>) ]N 
The direction cosines of the normal surface were defined in 
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equations (3.22), (3.23.3), (3.24) and the direction 

cosines of the ray at the front surface were defined in 
equation (3.26). The direction cosines of the skew ray at 
the rear surface, will therefore be: 



K * = n 12 


R z 


+ 


[cos (<j> 1 ) 


n 12 


cos ((f)) ] 


• n x 




L ' ’ n 12 


"y 


+ 


[ cos ( <f> ’ ) 


n 12 


cos ((f)) ] 


• n y \ 


(3.32) 


M' = n 12 


R z 


+ 


[cos (<f>' ) 


- n 12 


COS ((f)) ] 


• n Z ] 





3.2.6 Intercept at the Rear Surface 

The equation of the cone of the rear 



surface is 



2 2 2 
y + z = [ax - b ] 



(3.33) 



where a is the slope and is defined from figure 7, 



a tan(£j) _ x 



V 



R J + 1 R J 

X J + 1 X J 



(3.33.1) 



b is the intercept at the X - axis and is calculated from 
equation (3.33.1): 



b = X = 
v 



J + 1 



X J - R J 



\J + 1 



R J + 1 ~ R J 



(3.33.2) 



(XI, YI, ZI) is the point of intercept at the rear surface. 

The equation of the rear surface is therefore: 

YI Z + ZI 2 - ( aXI - b) = 0 (3.33.3) 

where a,b were defined in equations (3.33.1), (3.33.2). 
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\ 



\ 

\ 

\ 




Figure 6. Diagram used to calculate the direction cosines 

of skew rays at the rear surface. This diagram is 
based on Kingslake [16]. 
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A y 




Figure 7. Geometry of a truncated cone which represents a 
small portion of the rear surface of the lens. 
The cone angle is . 
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The direction cosines of the ray at the rear surface were 
defined previously (3.32). They are: 



K' = 



L ' = -- 



M' = 



ZI - zo 



(3.34) 



(XO, YO, ZO) is the point of intercept at the front surface 
(3.19.1, 3.19.2, 3.19.3). 

D = v/(XI - XO) 2 + ( YI - YO) 2 + (ZI - ZO) 2 (3.35) 

Solving for D from equations (3.33.3) and (3.34) one gets: 



D = 




(3.36) 



Where : 

a' = (L' ) 2 + (M' ) 2 - a 2 (K' ) 2 

b' = (L')YC + (M')ZO - (K ' ) a 2 (XO - b/a) 

c' = YO 2 + ZO 2 - a 2 XO(XO - 2b /a) - b 2 

Knowing D, the intercept point at the rear surface can be 
verified using equation (3.34): 



(3.36.1) 

(3.36.2) 

(3.36.3) 



XI = K'D + XO 
YI = L'D + YO 
ZI = M’D + ZO 



(3. 37.1) 

(3.37.2) 

(3. 37.3) 



3.2.7 Equation of the Normal to the Rear Surface 

The equation of the cone of the rear surface 
was defined previously (eq. 3.33): 

f(x.,y,z,) = ~(ax - b) 2 + y 2 + z 2 = 0 (3.33) 



21 



a,b are defined in equations (3.33.1), (3.33.2). 

7f = e + q §JL + e (3 38 1) 

X 5X Y 5Y Z 6 Z 1 J.JB.U 

Similar to what was done at the front surface: 



|| = -2a (ax - b) (3.38. 2) 

ff= 2y il= 2z (3.38.3) 

N- ># 2 + <H> 2 + 'ff * 2 ! 35 

~ r 2 - , i 2 , 2 2 

= 2 [ a ( ax - b ) + y + z ] 

= 2(ax - b) [ (1 + a 2 ) ^ ] (3.38.4) 



The direction cosines of the normal to the surface are as 
follows : 



_ (< t/ ao 

x " l?f| 



-a 



TPs 



(1 + a ) 



(3.39.1) 



n v = — 1 T-r 



(y" + 2 ") ' J d TT 7 )^ 

Z 



n z - 



(Y + z ) 2 (1 + a ) 






(3.39.2) 

(3.39.3) 



3.2.8 Refraction at the Rear Surface 

Similar to what was done at the front surface: 

R’ = K* (e x ) + L'(e y ) + M' (e z ) (3.40.1) 

^ A A A 

N = n x (e x ) + n Y (e y ) + n z (e 2 ) (3.40.2) 

Where the direction cosines of the refracted ray (R 1 ) and 

-V 

the normal (N 1 ) were defined previously (equations 3.32, 3.39). 
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The refracted angle at the rear surface: 

<f> T = cos tti (3.41) 

1 Ir'I . N'l 

Both the ray and the normal are normalised; therefore: 

= K'n x + L'n y + M'n z (3.41.1) 

4> t ' = sin ^t n 23 ♦ (3.41.2) 

3.2.9 Direction Cosines of Rear External Refracted 

gay 

As done previously (equation 3.32): 



K" = n 23 K' 


- [ COS (<+>-£*) 


- n 23 cos (^j) ] 


* n x > 




L" = n 23 L ' 


- [cos(4> I ') 


- n^cos ((fij) ] 




V (3.42) 


Z " = n 23 Z ' 


- [cos(^ I ') 


- n 23 cos ( j ) ] 


• n z 





3.2.10 Intersection with the Image Plane 



X lm " 3F = 


D 3 K" + XI ] 




Y Im = D 3 L " 


+ YI \ 


(3.43) 


Z lm = D 3 M ” 


ZI ) 





Solving equation (3.43), the intercept with the image plane 
can be calculated: 



D, = (BF - XI) /K" 



Im 


/ BF - 
' K 1 * - 


XI ' 


) L" 


+ YI 


(3.43.1) 


Im ~ 


/ BF - 

( k” 


XI 


) M" 


+ ZI ; 
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4. FORMULATION OF EQUATIONS FOR GRIN LENS DESIGN 



4.1 Solving the Integral for 9 



As mentioned in section 2: 

r 



9 e = 6 0 + E 



l 



r(rn(r) -E } 



(4.1) 



r^, r r are the radii to the front surface and to the rear 
surface, respectively. 



E = e . ^(r^) . r^ . sini^ 

n % ; ' A + B (%) 2 



(4.2) 

(4.3) 



Substituting: 

v 2 

v--(^) ; 

Equation (4.1) turns into the form: 



dV = _ 2r V 3 

dr r 0 



(4.1.1) 



9 = 9 n 

e 0 



E 

2r r 



v 



v. 



dV 



I 2 P 

AV + b - V -- 



rrp? 



(4.1.2) 



Therefore : 



9 = 

e 



-h 



_ 1 (2(E/r 0 ) V - A \ 

sin 



I'W 



<V r f> 



or: 



9 e = 



_ x ( 2E 2 /r 2 - A » 2E 4, /r f “ - A 

1 f -sin 



sm 



(A 2 + 4 BE 2 /r 0 2 ) ^ ■ 



,2 2 






(4.1.3) 
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The appropriate expression for the radius will therefore be: 



V 2 I E I 



(r, f ) 



v A + 7 A 2 + 4BE fc /r r J . sin{?29 Q + sin 



2 ^ 
/ v 

•o 



. -1 



2E;2 / r (f,r) - 



A 2 + 4BE Z /r Q 2 



(4.4) 



If the radius to the rear surface (r ) is calculated from a 
known radius to the front surface (r^) the minus sign is used 
in front of the angle 9 . This situation occurs with skew 
ray calculation. The plus sign is used when r^ is calculated 
from a known r . This situation occurs for lens design. 

4 . 2 Meridian Plane 

4.2.1 Solution for the First Ray 

(a) Define Point X^(l) , Y^( 1) 

As in the homogeneous case : 



R = 1.0 



Y^l) 
x x ( 1) 



R 

_R 

tana 



(4.5) 



(b) 

tan (9 



Define Angle to the Front Radius 
0 + 9 e ) X 1 (l) - X c 



(4.6) 



r f(D = r Q = 



sin ( 9 n ( 1) + 9 (1) ) 
u e 



(4.7) 



Refering to figure 8, 9 - is the angle to the rear radius, 9 
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is the angle between the rear and the front radii, and X c 

is the location of the center of symmetry. 

C c) Choose 0 from Initial Thickness 
e 

9 e ( 1 ) = T(l)/r f (l) 

9 0 ( 1) = (0 Q ( 1) + 8 e ( 1 ) ) - 0 e (l) 

(d) Solve for \p^ 

a + 1^(1) + u = tt/2 

1^(1) = tt/ 2 - a - u 

Knowing r^(l), the initial index of refraction (n 2 (r f )) 
can be calculated from eq . (4.3). Therefore: 

sind^(l)) = n 2 (r f )/n 3 . sin(I^'d)) 

I 1 ' (1) = sin -1 [n 3 /n 2 (r f ) . sin [1-^(1)]] 

Refering to figure 8, \p^ can be calculated. 

9 e (1) + 0 O (1) = a + 5 

•*•6 = 0(1) + 9 n ( 1) - a 

e u 

5 + ip f (l) - I 1 ’ (1) = tt/2 
Tp f d) = TT/2 + I 1 ' (1) - 6 

(e) Solve for E 

E = e . n (r-) . r c . sintd 

r f f f 



(4.8) 



(4.9) 



(4.10) 



(4.11) 



(4.12) 



(4.2) 
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Normal 
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Figure 8. Diagram degining angles and coordinates for a ray transmitted by the 
GRIN lens. For the case illustrated B is negative. 



Where : 



e = ±1 



ip < ir/2 ■* £ = +1 

ip > tt/ 2 e = -1 



(4.2.1) 



(f) Solve for r f ,^ r 



Knowing the radius to the front surface 



(r f ) and the angle between the two radii (6 e ) , the radius to 
the rear surface (r^) can be calculated using equation (4.4). 
Using equation (4.2), \p can be calculated: 



( g) Solve for 8 ( 1) 

The focal length, which is an input, is known 
BF = F/R 

The coordinates at the rear surface are: 




(4.13) 



X 2 (l) = X c + r r cos(e Q (l)) / 
Using equation (4.13), 8(1) can be calculated: 



8(1) tan bf~- x~TI5" 



-i r y 2 (i > 



(4.14) 



(h) Solve for the Refractive Angles at the 



Front Surface 



As in the homogeneous case: 
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n (J) + i 2 ' (J) = 


B (J) 


(4.15) 


sin[I 2 (J)J = n 2 /n 2 (r r ) . 


s in [ 1 2 ' (J) ] 


(4.16) 


ip (J) + I,(J) + n ( J) 


+ e n (j) = tt 


(4.17) 



Solve for n (J) in a manner similar to the homogeneous case: 



Where : 



n ( J) 



sin 



-1 



s/ 1 + a ' 



A 



cosS' - n 2 /n 2 (r r ) . cos0(J) 
sin<S' - n 3 /n 2 (r r ) . sinS (J) 



(4.18) 



(4 . 18. 1) 



Where : 

<S' = it - e Q (j) - ipg (j) 



(4.18.2) 



Knowing n(J) , I 2 ' (J) can be solved from eq. (4.15) and then, 

I 2 (J) can be solved from equation (4.16). 

4.2.2 Solution for Next Rays 

4 . 2 . 2 . 1 Coordinates at the Rear Surface 
Solve for the subsequent rays the 

same way as in the homogeneous case (equations 3.10, 3.10.1). 

4 . 2 . 2 . 2 Calculate Radius and Angle .to the 
Rear Surface 

r r (J) = y ( X 2 (J) - X c ) 2 + Y 2 (J) 2 (4.19) 



6 0 (J) = tan _1 [Y 2 (J)/(X 2 (J) - X c ) ] 



(4.20) 



29 



4. 2. 2. 3 Solve Front Surface Parameters 



Refer to Figure 8 for a definition 

of geometrical variables. 

(a) Equations 

n(J) + I 2 ' (J) = B(J) (4.15) 

sin[I 2 (J)] =n 3 / n 2 ^ r r ^ * (4.16) 

lp r (j) + I 2 (J) + T) ( J) + 0 0 ( J) = TT (4.17) 

E = e . n 2 (r r ) . r r . sin[^ r (J)] (4.2.2) 

E = e . n 2 (r f ) . r f . sin[^ f (J)] (4.2) 

0 (J) + 0 n ( J) = ct + <5 (4.21) 

<S + 4> f (J) - I 1 ‘ (J) = tt/2 (4.22) 

r f = £(r r' 0 e (J)) (4,4) 

s in [ I ^ ( J ) ] = n 2 (r f )/n 1 . sinflj^fJ)] (4.23) 

a + I-^(J) + u = tt/2 (4.24) 

(b) Solution by Iteration 

An explicit solution is not 

possible. An iterative solution is necessary. 

(bl) Assume 6 (J) . 

e 

(b2) Solve for coordinates 

of the front surface: 

- Calculate I^(J) from eq. (4.24) 

- Calculate 6 from eq. (4.21) 

- Calculate front surface coordinates and radius: 
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(4.25) 



Y X CJ) = tan(9 0 (J) + 0 e (J)) . (X^J) - X c ) 

Y 1 (J) = tana . X 1 (J) 



The result is: 



X-^J) 



Y 1 (J) 



X c . tan[0 Q (J) + 0 e (J) ] 
tan [ 0 q ( j 5 + _ G ( J ) ] - tana 



tana . X^(J) 



r f = Y 1 (J) /sin[e Q (J) + 9 e ( J ) ] 



X c . tana 



(4.26) 



(4.27) 

(4.26) 

(4.28) 



sinleT q Tj) + 9 e ( J FI - tana I cos fQg ( J) + 6 (j) ] 

(4.28.1) 



(b3) Knowing r . , 1^' (J) can 
be solved using equation (4.23) . 

- Calculate ip ^ from eq. (4.22) 

- Calculate E from eq . (4.2) 

(b4) Calculate second expression 
for front radius (r^) from eq . (4.4). 

(b5) The iteration routine: 

The iteration search 

for correct 8 g uses the Newton- Raphson routine. 



The function: 



F = r , r - r, = 0 



(4.29) 



r^ is the expression from the formula (eq. 4.4) where x c is 

the geometric expression (eq. 2.28). 

The derivative is found with resDect to 9 : 

e 
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F' 



(4.30) 



d(r f f ) d(r f ) 

= dToj dTej 
e e 



d(r ff> 

dTO" 



= + 



^ r ff \ . cos[+20 e + sin ^ ( B 2 ) ] 

2E 2 



(4.30.1) 



Where : 

B x =J A 2 + 4B(E/r r ) 2 (4.30.2) 

B 2 = [2(E/r r ) 2 - A ] /B x (4.30.3) 

The lower signs are used when the front radius (r f ) is calculated 
from a known rear radius (r r ) (this case occurs at J > 1) , while 
the upper signs are used when r r is calculated from r f (for J=l). 



d(r f ) -X c . tana . {cos[9q(J) + 9 e (J)] + tana . sin[0 Q (J) + 0 g (J 
jsin [ 9g ( J) + © e ( J) ] - tana . cos[0 Q (J) + 6 e (J)]| 2 



(4.30.4) 



(c) Correct Solution 

After knowing the correct value 
for 0 g (J), the procedure used in (b) can be followed resulting 
in the correct values for: 



I 1 (J), 5, X 1 (J), Y X (J), r f (J) , I 1 '(J), r{> (J) , E. Using 
equation (4.2.2) $ can be calculated. Equations (4.15), 

(4. 16) , (4. 17) can be solved for n(J), I->'(J), B(J) as shown 
previously for the J=1 case (equations 4.13, 4.15, 4.16). 
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4.3.1 Front Surface 



As in the homogeneous sace , the intercept 
point at the front surface (XO, YO , ZO) can be found (equations 
3.19.1, 3.19.2, 3.19.3). 

Also, the direction cosines of the normal to the front 
surface can be calculated (equations 3.22,3.23). The direction 
cosines of the ray can be found as well using equation (3.26). 

4.3.2 Front Radius 

Knowing the intercept point at the front 
surface, the front radius can be calculated: 

r f (J) = [(XO - X c ) 2 + (YO) 2 + (ZO) 2 ] 55 (4.31) 

Knowing r^.(J) , the initial index of refraction [n 7 (r^.)3 can 
be calculated using equation (4.3). 

4.3.3 Refraction at the Front Surface 

As for the homogeneous case (section 3.2.4), 
the refracted angle at the front surface can be calculated 
using equations (3.30), (3.31). 

4.3.4 Initial Direction Cosines for the Internal 

Refracted Ray 

The initial direction cosines of the internal 
refracted ray (K',L , ,M' ) are calculated the same way as for 
the homogeneous case (section 3.2.5, equation 3.32). The 
vector form of the initial direction of the ray is therefore: 

/s /v /v 

R = K' . e„ + L ' . e v + M' . e„ (4.32) 

A i ci 
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4.3.5 Front Radius Direction Cosines 



The vector form of the front radius is 
defined as follows: 

_y A A A 

r f = r fx • e x + r fY * e Y + r fz * e Z 
Where the direction cosines are: 

r fX = (XO - X c )/r f (J) 
r f y - Wr f (J) 
r fz = ZO/r f (J) 



(4.33) 



(4.33,1) 



4.3.6 Normal to Plane of the Curved Ray at (XO,YO,ZO) 
The vector cross product of r^ and R is normal 
to the plane defined by the two vectors. 



NP f = r f X R = 



'X 



*fX 
K' L' 



'X 



r fY r fZ 



M' 



(4.34) 



Where : 



NP^ ' e X + ^ P fY * e Y + ^ P fZ * e Z 



(4.34. 1) 



and 



NP fx - 


r fY 


. M- - r fz . 


L ' \ 




NP fY = 


r fz 


. K< - r fx . 


M' > 


(4.34.2) 


NP fZ = 


r fx 


. L' - r fy . 


K' / 




NP f =V 




+ NP fy + NP 


2 

fZ 


(4. 34.3) 



The normalized direction cosines of the normal to the plane 
of the curved ray at (_XO,YO,ZO) are therefore: 

NPfx/NPf ; NP fy /NP f ; NP fz /NP f 
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(4.34.4) 



(4.35) 



4.3.7 Calculating the Angle Between the Ray and 
the Front Radius 

r f . R = jr f | . R| . cos(i|j f ) 



Both the ray and the unit vector for the front surface radius 
are normalized, i.e. 



hi 

Hi 


, 2 
* (r fx 


+ r 2 + r 2 )^ 

r fY fZ ; 


R 


CM 

W 

II 


+ L' 2 + M' 2 ) ^ = 



(4. 35. 1) 



There fore : 



ipf = cos [r^ . R] 



= cos -1 [r fx . K' + r fY . L' + r fZ . M' ] 



(4.35.2) 



4.3.8 Calculating constant E 

Knowing constant E can be calculated 

using equation (4.2) . 

4.3.9 Intercept Point of the Curved Ray with Rear 
Surface (Point XI,YI,ZI) 

4 . 3 . 9 . 1 Equations 

As previously, 0 e is the angle between 
the front radius and the rear radius. 



r rX * r fX + r rY * r fY + r rZ * r fZ = cos £ 9 e (J)] (4.36) 

The direction cosines of the front radius were calculated in 
equation (4.33). Similiarly, the direction cosines of the 
rear radius are: 
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r rX ' (XI - V r r (J) ) 

r rY = YI/r r (J) > (4.36.1) 

r rZ = ZI/r r (J) ' 

Therefore, equation (4.36) can have the form: 

(XI - X c ) r fx + (YI) r fY + (ZI)r fZ = r r (J) . ccs[9 e (J)] (4.36.2) 
The rear radius is defined as: 

(XI - X c ) 2 + (YI) 2 + (ZI) 2 = [r r (J)J 2 (4.37) 

The rear surface was defined previously. Consequently: 

- (aXI - b) 2 + (YI) 2 + (ZI) 2 = 0 (3.33) 

a,b were defined by equations 3.33.1 and 3.33.2. 

Also, from geometry: 

(XI - X c ) = r r (J) . cos [6 q ( J) ] (4.38) 

Another equation arrives from the connection between the two 
radii through eq. (4.4) . 

4 . 3 . 9 . 2 Solution 

a. Assume 6 (J) . 

e 

b. Calculate r r (J) from the formula 
(equation 4.4). 

c. Calculate XI from equation 
(4.36.2) : 

XI = r r (J) . cos [ 9 g ( J) ] + X c (4.36.3) 

d. Calculate radius from equations 
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(4.37) and (4.36.3) : 



Radius = R = [ (YI) 2 + (ZD 2 ]* 5 

= r r (J) [1 - cos 2 [0 Q (J) ] ] 

= r r (J) . sin [ e Q ( J) ] (4.37.1) 

e. Calculate parameters 

a and b from equations (3.33.1) and (3.33.2) as in the 
homogeneous case. 

f. Substituting into equation 
(3.33) : 

XI = (R + b)/a = (r r (J) . sin [0 Q ( J) ] + b }/a (4.38) 

where R is the radius defined in equation (4.37.1). 

g. Solving equations (4.36.3) , 
(4.38) for 9 q(J) one gets: 

/(B1L 2 - (AH (CM 1 ,4.39) 

A 

Where : 

A' = [r r (J) ] 2 (1 + a 2 ) 

B ' = [r ( J) ] (b - a . x c ) 

C' = (b - a . X c ) 2 - a 2 . [r r (J) ] 2 

h. Knowing 0 q(J) , it is pos- 
sible to calculate new value for 0 (J) and iterate through steps 
(a) - ( g) until the difference in the values of 9 q(J) is 



(4.40) 



0 Q ( J) = sin 



-1 



3 ' 

'IT 
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small enough. 



i. When an accurate value of 

0q ( J) is received, XI is calculated from equation (4.36.3). 

YI and ZI are calculated using equations (4.36) and (4.37.1). 
Consequently : 



YI 

ZI 



. B" , v(B" ) 2 - (A") (C") 
A ,r A" 



= N/{r r (J) . sin [ 0 Q (J) ]} 2 



(YI} 2 



Where : 



A" 

B" 

C" 



1 + [r fy ] 2 /[r fz ] 2 

(r r ( J) . cos [0 g ( J) ] 
(J) 

. cos [ 0 g ( J) J 

fZ 




(XI - X c )r fx )r fY /r fz 
<XI - V^| 2 - j r r <J > 



When r^ z = 0 , 



the solution for YI and ZI is different: 



(4.41) 



sin [ 0 q (J) ] 
(4.42) 



YI = (r r (J) . cos[0 e (J)] - (XI - X c )r fx }/r fY (4.43) 

While ZI is calculated by equation (4.41). 

4.3.10 Direction Cosines of the Ray at the Rear 
Surface 

> a. Similarly to what was done at the front 

surface (equation 4.35): 

r r . R = jr r j . R . cos(^) (4.44) 

Both the ray and the rear radius are normalized; therefore: 

r rX • K r ' + r rY * L r' + r rZ * M r ' = cos ^ r (4.44.1) 
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b. The ray remains within the same plane 
while transfering from the front surface to the rear surface, 
i .e . : 



NP . K ' + NP W . L ' + NP,„ . M ' = cos ( ir/2 ) =0 (4 

fX r fY r fZ r 



normalized, i.e. 



c. The direction cosines of the ray are 



(K r ') 2 + (L r M 2 + (M r ') 2 = 1 



U 



Where : 



d. Solution: 

(dl) NP fx f 0, NP fy 1 0, NP fz ft 0 



(L V 


B" . y/(B”) 


2 


~ 2A !r ~ 


2 A 


(M' r ) 


- — - — (L' ) 

3 B L r’ 




(K' r ) 


= [-NP fY . (L' 


J - 


A = NP fx 


. r rY - NP fy . 


r rX 


B = Np fx 


. r rZ - NP fz . 


r rX 



. " \ / n >' 



fZ 



U 



(M' r ) ]/NP fx 



C = NP fx . cos(^ r ) 



A' = 


( NP fY 




Up fx 


B' = 


/ NP fZ 


1 NP fX 




2 (NP 


C = 





) 



+ 1 



2 

) + 1 



fY ; 



(NP fx ) 



fZ' 

T~ 



(< 



. 45 ) 



.46) 



1.47) 



.48) 



39 



A" 



B" 

C" 



(A*) + (B') (|) 2 
2 (A) (C) (B 1 ) 

(Bp 

(BM (|) 2 - 1 



(CM (A) 

(§T“ 

(O’) (C) 



(4.48) 



In (4.47) the (+) sign is used when r Y < 0 while the (-) 
sign is used when r rY > 0 . 

(d2) In the above case, when B = 0 (look for 
definition of B in equation (4.48)) the solution for direction 
cosines of the ray is different. 



(L' r > 



c 

A 



(M' r ) 



(C) (L’ r ) 



/ r (C') (L' r ) , 2 



\ / 



“(3') [(A*) (L' r ) - 1] 



B 



(K’ r ) = [-NP fy . (L' r ) - NP fz . (M’ r )]/NP fx 



(4.49) 

All the symbols used here were defined in equation (4.48) and 
the ± sign is used as before. 

(d3) NP fx = 0, NP fy = 0, NP fz ^0, r rY ^ 0 

(M' r ) = 0 (4.50) 
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R ** 

(K' ) = — + 

^ r A" - 



(B") - (A") (C") 

A n 



(L V = < cos( V - r rX ' K 'r ! 



The (+) sign is used here. 
While in this case: 

A" = 1 + (In) 2 
' r rY 



(4.50) 



B" = cos(^ ) — 

< r rY> 



(4.51) 



C" = 



cos (to) 



rY 



- 1 



(d4) In the above case, 



soluticn is: 



(M' r ) = 0 

cos (to ) 

(K 1 ) = — 

r* r 

rX 



(L' r ) = ±Jl ~ (K ’ r ) 2 



Again, the (+) sign is used here. 

(d5) NP fx = 0, NP fy jt 0, NP fz =0, r r2 ? 0 



The solution is siirilar to that in case (d3) . 



when r ^ = 0 the 



(4.52) 
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(L V 



(K' r ) 



(M* ) 
r 



0 



B" x V(B") 2 - (A") (C") 

A ir 1 



““ [ COS ( 'p ) 

rZ 



rX 



K' 3 
r 



(4.53) 



(A"), (B"), (C") are the same as defined in (4.51) except that 

r rV is replaced everywhere by r rZ . The (+) sign is used here. 

(d6) In the above case, when r r2 = 0 the 
solution is similar to that in case (d4) : 



(L' r ) = 0 



cos (ij; ) 



(K' r ) = - 



rX 



(M' r ) = ± s/l - (K' r ) 2 



(4.54) 



Also here, the (+) sign is used. 



(d7) NP fx = 0, NP fy ji 0, NP f2 jt 0, r rX = 0 



(L ' r » - - 



cos ( ) 



NP 



fY 



’rY NP * rZ 



NP 



(M'_ ) = - 



fY 



(L ' ) 



(4.55) 



fZ 



(K’ ) = ±V 1 - (L' ) 2 - (M‘ ) 2 
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The ( + ) sign is used here. 

(d8) As in (d7) , but: (r 



NP 



fY 



rY NP 



r v jt 0 
rX r 



fZ 



r rZ> 



= 0 ; 



(K' ) 
r 



cos (\p ) 
r 



rX 

r 



(L’ r ) = ± 



1 - (K' r ) 2 

77CW 



(4.56) 



NP 



fY 



(M V =- NP^ (L V 



(d9) As in (d7) , but: r rX ^ 0 



, _ (B") . V(B") 2 ~ 4 ( A ") ( C") 

(L r 2TPT 2 (A 11 *) 



NP 



(M‘ ) = 
r 



fY 



NP 



(L* ) 



fZ 



NP 



(K' ) = [cos (to) - (r y 

rX 



fY 



NP f g r rZ } (L ' r } ^ 



(4.57) 



As in case (dl), the (+) sign is used when 
the (-) sign is used when: r rX > 0. 



r rY < ® while 
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5. FLOW CHART OF PROGRAM LENS 



(Homogeneous Case) 

Note: Subroutines SHAPET, SKEWT compute the shape of lens 

and the skew rays as described in the previous sections (3,4). 
However, subroutines SHAPEC, SKEWC compute the shape and the 
skew rays of a different design of the lens. In this design, 
the fixed cone is the rear surface instead of the front 
surface (compare figure 9 to figure 2.). The equations for 
this different design are not included in the current report.. 
The equations are part of Capt. Carr’s thesis [15] and will 
be described there. 
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Rear Surface 




CD 



3 

(f) 



C 

o 
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Figure 9 . Illustration of the lens which has a circular cone on the inside 
surface . 



P roar am LENS 



Main 
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47 




n 1 ,n 2 ,n 3 



T (1) = T tl)/B 



IT, RAD 
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In program GRIN, which uses the same subroutine 
data for XC, DIV is given. 



here 



SHAPET 








Cor 

ay 

dx 


npute 

rr (j) 

T 


! 


Co 

dv 

dx 


mpute 

N (J) 



^ Return 
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POINTS 



Compute 

x^(J) 



Compute 

v n ( J) 



Compute 
x 1 (J) 



Compute 
y : (J) 

Compute 



T ( J) 



C 



Return 
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CALETA 



i 



) 



Compute 

n(J> 




54 



SHAPEC 



Compute 

8 ( 1 ) 
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0 




56 




57 







(*) 



In program GRIN, which uses the same subroutine, here, 
also >:c is shifted. 




Ntncty(l) = 1 





Compute 
Ceni ,Torbi, 
Boti 



© 
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60 




61 
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© 



62 




63 



RADI 




© 




65 



0 




66 
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© 




( 



Return 






SKEWT 



Compute 
xO , yO , zo 



N —sin (a) 
x 





@ 



0 




0 
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© 




I 
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Compute 
K", L", M" 



Compute 
Dl, D3 



Compute 

Opl 



Compute 
y iM (IZ) ' Z IM (I 



© 
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© 



PRINT 




©ET - 

Ntncty (Iz) =1. 



Define 
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© 



ICHECK=0 
Ntncty (Iz) =0. 




© 



ICHEC 
Ntncty 1 


:k=o 

(Iz) =0. 






Flag(Iz ) —2 







PRINT 




© 



© 



ICKEC 

ntncty 


:K=0 
(Iz) =0 . 






Flag(Iz) =3 







PRINT 
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© 





7 G> 
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Rays=Iz-l 



Drays=Rays-Iz^ 




Coord = 1 

E1 J )= ~ y l° 
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Compute 
ElyO. E1 z2 



PRINT 




Coord=Coord 
+ 1 

El () =E1 +ELL 
z zz 




YES 

© 
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CoorO =Coord 
Coord =Coord 



© 
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Define 
Nuinfcer ,Numb 




| YES 






© 




82 





83 



© 




84 





35 



PLOTT 





36 
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6. FLOW CHART OF PROGRAM GRIN 



Note: Program GRIN is identical to program LENS (Section 5) 

except for subroutines: 

SHAPET, SKEWT, SHAPEC, SKEWC. 

The flow chart for the first two is listed below. Subroutines 

SHAPEC, SKEWC are of a different design (fixed cone a for 

P 

the rear surface). They are part of Captain Carr's thesis [ 15 ] , 
and they will be described there. 
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© 
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91 




92 




c 



Re turn 



) 



9 3 





94 




95 




© 



96 






97 



0 




98 



0 




99 




100 



© 





Return 



INDEX 







Compute 

N 2 







( Return ) 
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Dif f= 

10 (NEW) -0 
e e 

(OLD) I 
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DIRECT ) 










Return 



© 






10 7 




© 




10 8 




c 



Return 



J 
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0 







Return 



7. COMPUTER PROGRAMS USER'S GUIDE 




case 



a , Grid, Ellnum 

fr 

X c , Div 

Shape, Squray, Ellips 

7 . 2 Options 



Rear surface is a fixed cone 

CaS6 

(Capt. Carr's thesis, see ref. [15]) 

I =2 Front surface is a fixed cone 

case 

(This report) . 

7 . 3 Execution Commands 

7.3.1 Opening Commands 

LAXXXX (Write account no. instead of XXXX) . 
Password 



GLOBAL TXTLIB FORTMOD2 MOD2EEH 

7.3.2 Coinpilation 

FORTGI LENS (or GRIN) 

7.3.3 Run on Terminal (WORKG^) 

FILEDEF 05 DISK DATA L 

FILEDEF 06 DISK OUT D (RECFM F 3 LOCK 132 PERM 
EXEC RUN LENS (or GRIN) 

XEDIT OUT D 
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7.3.4 



Hard Copy 



FILE 

PRINT OUT D 

7.3.5 Get a Plot 
LENSCOM PRINTER 

7.3.6 Change in Plot Routine 
XEDIT LENSPLOT FORTRAN 

7.3.7 Transfer Files to Another User 

(a) The Transferer 

CP SPOOL PUNCH CONT TO XXXXP (2) 

PUNCH LENS FORTRAN (3 ) 

(PUNCH GRIN FORTRAN (3) ) 

(PUNCH DATA L (3j) ) 

CP SPOOL PUNCH CLOSE NOCCNT 
CP SPOOL PUNCH OFF 

(b) The Re c eiver 
READCARDA *A * 

(c) Notes 

(1) WORKG is the name of the file of the 
commands in this section. 

(2) Account No. of the receiver. 

(3) Optional. 
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8. COMPUTER PROGRAMS LIST OF SYMBOLS AND SUBROUTINES 



8.1 


Subroutines 


DATTA 


Initial data 


SHAPET 


Computes shape of lens in I =2 

C cLS 0 


SHAPEC 


Computes shape of lens in I =1 

case 


INIT 


Computes behavior of first ray. 


POINTS 


Computes coordinates of rear and front surfaces. 


FRONT 


Computes I^(J) , 1^' (J) ,ty,E at front surface. 


REAR 


Computes n ( J) ,^ 2 (J) >*2 ' at rear surface. 


CALTC 


Computes thickness at each point for SHAPEC. 


NOSE 


Computes extrapolation of last points in I =1, 

CcLSG 




or coordinates shift in I =2. 

case 


ZOYO 


Computes Zq ' (IZ) , y Q ' (IZ) for skew rays calculation. 


RADI 


Computes intersection of the ray with the designed 



surface . 



SKEWT 


Skew rays calculation in I =2. 

esse 


SKEWC 


Skew rays calculation in I =1. 

esse 


INTENS 


Calculates transmitted intensity of the rays. 


IMAGE 


Calculates image plane spot diagram. 


SUMM 


Calculates spot diagram energy density. 


PLOTT 


Save output for plot routine. 


INDEX 


Computes index of refraction in the lens (NO . 


RADN2 


Computes radius to rear (or front) surface from 



known radius to front (or rear) surface and e g . 



ITER 


Iteration routine to find correct 6 g in SHAPET. 


DIRECT 


Computes direction cosines of the ray at the rear surface 
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8.2 


Program Symbols 
8.2.1 Data Symbols 


F 


Focal length . 


a 


The fixed surface cone half angle. 


u 


The incident ray offset angle. 


I 


Number of iterations. 


R 


Maximum radius of fixed surface cone. 


T ( 1) 


Thickness of lens at the edge. 


n 1 ,n 2 ,n 3 


Indices of refraction. 


a 

P 

Grid 


Inclination of object plane. 
Size of grid increment. 


Ellnum 


Number of ellipses. 


Ell 


z - coordinate increment for ellipse plot. 


X C 


Location of center symmetry in GRIN. 


Div 


The maximum change in the index of refraction 

<v 

8.2.2 Geometrical Symbols 



x 1 (J),y^(J) Coordinates on front surface. 

X 2 (J),y 2 (J) Coordinates on rear surface. 



r f (J) 

-L- 


Radius to front surface . 


r r (J) 


Radius to rear surface • 


XO , YO , ZO 


Intercept of the ray with the front surface • 


XI , YI , ZI 


Intercept of the ray with the rear surface • 


Y_ ,Z T 
Im' Im 


Intersection with the image plane • 
8.2.3 Direction Cosines Symbols 


N ,N ,N 
x' y ' 2 


Direction cosines of the front surface normal 
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N t ,N t# N _ 
xl yl zl 

R / R , R 
x y z 

K' ,L 1 ,M' 



r £x' r fy ' r fz 

NP f x,NPfy,NP fz 

r rx' r ry ' r rz 
K" ,L" 



8 . 2.4 



Direction cosines of the rear surface normal. 
Direction of cosines of the ray. 

Direction cosines of the internal refracted 
ray. 

Direction cosines of the front radius. 

Direction cosines of the normal to plane of 
the ray at front surface. 

Direction cosines of the rear radius. 

Direction cosines of the rear surface external 
refracted ray. 

Angle Symbols 

The angle symbols are defined in figures 3 and 8. 
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9. PROGRAM LENS 



LISTING 
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FILE: GRIN FORTRAN A NAVAL POSTGRADUATE SCHOOL PAGE 019 



11. PRELIMINARY RESULTS 



Table 1. Comparison of RMS Spot size for a homogeneous lens 
with the RMS Spot size for a GRIN lens. 



X ' S X s piV 

a p 


0 


+0.05 


in 

o 

• 

o 
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0 . 


0.0057336 


0.0069021 


0.005286 


0.1 


0 .027424 


0.027646 


0.027327 


0.2 


0.06277 


0 .06288 


0.06304 


0.3 


0. 113431 


0.11216 


0.113245 


0.4 


0.18265 


0.18287 


0.18277 



Note: DIV is the fractional variation of refractive index, 

cip is the angle of incoming parallel rays relative to lens 
axis; otp has dimensions of radians. The center of symmetry 
X c was at -0.0 7 measured from the apex of the outer surface. 
The minus sign indicates that X r is external to lens. 
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RESULTS FROM LENS 



(ICASE = 1) 

Input: a = 0.0, 0.2, 0.4 
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GLM LENS SHAPE 
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0000000 T = 0.0500000 INDICES OF REFRACTION ARE: N1 = l .0000000 N2 * 1.5000000 N3 = 1.0000000 
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Gl M LENS SHAPE 
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I MAGE PLANE 
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SUMMARY HF PFSULTS 



S LINO :^'IVDS 



i 



i 

o 



i 

o 



c 



o 

w 



o 

o 



as 

o 



• *r 



***♦«■*♦ *********f** ♦ + ' 



Q 4- 

^ fr 



C < • 

• -f • 

c + • 

£ 4 • 

i • 



• • 

• • 

• • 
• • 



I 

**- 

*• 

. 4-0 
♦ IV) 



♦ 

♦ 

♦ 

♦ 

♦ 



• • 

• » 



♦ 

* 

♦ 

* 

♦ 

* 

♦ I 
♦o 
♦• 



♦ 

♦ 



♦ 

*o 

tc 

♦'0 

♦ 



4* O' 



♦ 

♦ 

♦ 

♦ 

•* 

♦ 

♦ 



c * 



*• 
+ rvj 
♦O 
■ * 



l 

o 



c 



o 

U) 



UJ 

U' 



o 

-J 

o 



■e- 

u 



CD 

o 



170 



SCAIF: "«"• = 0.S18F-01 UNITS 




* 

♦ 

*• 



o * 
• ♦ 

o * 



***#«***«*>*‘ 



k * * * * ♦ ***•♦**♦#•*♦ * #•*•***- 



. *-*•*•■• . 



I I I 



rvj 


r\j 


►— 


• 


• 


• 


u> 


o 


vr 


c 


-*< 


VJl 



I 



o 



I 

o 



\r 

N> 



o 



o 



171 



/v /> /’X 



~\ -\ 



A -X A 






JfJVld j t * V <* 1 



Siiun I 0-JO?1 “0 = IV 



o 


►— 


N» 


u> 




• 


• 


• 


• 


• 


o 


ro 


■f* 


o 


CD 




O 


o 


o 


O 




o 


o 


o 


a 




T1 


m 


Tl 

I 


•Tl 

1 




1 

o 


o 


1 

o 


1 

o 



I 

o 



O ♦ **-****.* + .♦. #*♦*■**** 



► ♦***- 4 *-*#*-* + *■****+**♦ + #*+***• 



*• • 



a 

_n 



I 

c. 



0 * 

— *. 

1 * 



# 

♦ 

'V * 



c * 
o * 



rr, » 

I > 



vC * 

w <► 



rs, * 
o * 
c_ 

m * 



■1* *“ + ♦****♦*<■ 


♦ * ♦ ****♦■*♦ * -ti- 


#♦■** **«* 


#■*♦■*-*•* + 


o *~ 


ro 


w 




• • 


• 


• 


• 


o N> 


X- 


O' 


CD 


o 


o 


o 


O 


c 


o 


c 


o 


m 


m 

i 


m 


m 


c 


o 


o 


O 


r— 


i— 







172 
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RESULTS FROM GRIN 



(ICASE = 2) 



Input : 



X c = -0.07, DIV = +0.05 

a = 0 . 0 , 0.1 
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1. Checkup for the homogeneous case. 
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R ESULTS FROM GRIN 

2. Input: X c = -0.07, DIV = -0.05, a p = 0.0 
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12. SUMMARY OF THE REPORT 



Lens design and ray tracing to determine optical 
aberrations requires efficient computer cedes. A clearly 
structured computer code with efficient algorithms has been 
developed to design a lens incorporating gradient refractive 
index, GRIN. One surface of the lens, either the inside or 
outside surface, is a circular cone. 

The computer code includes lens design and calculation 
of optical aberration for a lens using a homogeneous material, 
i.e., for n equal to a constant. The homogeneous case is 
included for several reasons. The simpler case with n equal 
to a constant provides insight to the more complex case of 
variable n. Extensive computer codes require a check on the 
validity of the output. The homogeneous case provides a check. 
When the factor B is equal to zero in equation [2.9], the 
results for the two computer programs with and without GRIN 
should be identical. Needless to say, the results are identical 
for the programs reported here. 

The case of a lens with GRIN and with a circular cone for 
the outside surface has been programmed. The companion case 
with the circular cone on the inside is being developed by 
Carr [15]. Carr's program can be readily incorporated into 
this computer code. 

The program has been exercised for the case of 5% 
variation in refractive index for both negative and positive 
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values for B. The aim of the preliminary calculations 
was to verify operation of the computer program. More extensive 
calculations are required in order to understand the influence 
of various variables. As more experience is gained, the 
program should be combined with an optimization code to 
permit thorough and efficient lens design. 

The selection of a 51 variation in refractive index 
was based on current state-of-the-art in the production of 
GRIN materials. A wider variation in refractive index may 
be needed. 

The program is restricted currently to values of X c 
satisfying the constraint X c < 0 . Additional work is needed to 
permit designs using positive values of X c « When X c is 
positive, the angle 0 at some point along the ray interior to 
the lens has a critical value at which ar/d6 changes sign. 

The critical value causes computational difficulties. 
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